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Abstract 

A technique for robust identification of nonlinear dynamic systems is developed and 
illustrated using both digital simulations and analog experiments. The technique is based 
on the Minimum Model Error optimal estimation approach. A detailed literature review 
is included in which fundamental differences between the current approach and previous 
work is described. The most significant feature of the current work is the ability to 
identify nonlinear dynamic systems without prior assumptions regarding the form of the 
nonlinearities, in contrast to existing nonlinear identification approaches which usually 
require detailed assumptions of the nonlinearities. The example illustrations indicate that 
the method is robust with respect to prior ignorance of the model, and with respect to 
measurement noise, measurement frequency, and measurement record length. 

Introduction 

The widespread existence of nonlinear behavior in many dynamic systems is well- 
documented, e.g, Thompson and Stewart [1]; Nayfeh and Mook [2]. In particular, 
virtually every problem associated with orbit estimation, flight trajectory estimation, 
spacecraft dynamics, etc., is known to exhibit nonlinear behavior. Many excellent 
methods for analyzing nonlinear system models have been developed. However, a key 
practical link is often overlooked, namely: How does one obtain an accurate mathematical 
model for the dynamics of a particular complicated nonlinear system? 

Accurate dynamic models are necessary for analysis, filter design, and/or control 
system design. For example, most filter design assumes white process noise, yet many 
nonlinear effects are inherently non-zero mean; e.g., quadratic nonlinearities are always 
positive. In order to obtain a model with truly zero mean process noise for filter design 
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purposes, all of the quadratic terms (and many other nonlinearities) must be well modeled. 
However, the complexity of many real systems precludes the possibility of accurately 
constructing a dynamic model purely from analysis using the laws of physics. 

Identification is the process of developing an accurate mathematical model for a 
system, given a set of output measurements. Much work has been done on identification 
of linear systems, resulting in a number of efficient algorithms. The accuracy and ease of 
application of these algorithms has given linear identification an enormous popularity. It 
is, therefore, a common practice to use linearized models to describe nonlinear systems. 
However, linearization does not work in every application, and even when it does provide 
a reasonable approximation, the approximation is normally limited to a small region about 
the operating point of linearization. Consequently, there is a real need for nonlinear 
identification algorithms. If nonlinearities are a predominant part of a system s behavior, 
using a linear model to describe such a system leads to inconsistencies ranging from 
inaccurate numerical results to misrepresentation of the system’s qualitative behavior. 
Since nonlinearities are seldomly easily characterized, identification techniques may prove 
beneficial in developing accurate mathematical representations of nonlinear systems. 

Numerous methods for the identification of nonlinear systems have been developed 
in the past two decades (Natke, Juang and Gawronski [3]). Most methods fall into one 
of the following categories: 

1. describing the nonlinear system using a linear model 

2. the direct equation approach 

3. representing the nonlinear system in a series expansion, and obtaining the respective 
coefficients either by using a regression estimation technique, by minimizing a cost 
functional, by using correlation techniques, or by some other approach 

4. obtaining a graphical representation of the nonlinear term(s), then finding an analytical 
model for the nonlinearity 

With such diversity of nonlinear identification techniques, the choice of an algorithm 
may be based on criteria such as: iterations required, robustness in the presence of mea- 
surement noise, number of measurements needed, robustness with respect to knowledge 
of the inital conditions, and robustness with respect to initial assumptions regarding the 
form of the nonlinearity, depending on the needs of the particular application. 

Among the methods which linearize the nonlinear system are those presented by 
Jedner and Unbehauen [4] and Ibanez [5]. Jedner and Unbehauen [4] represent a nonlinear 
system, which may often function at a number of operating points, by an equivalent 
number of linear submodels. It is assumed that the system operates at only a few points. 
Although the model is good for controller design, the point at which the system is 
operating must be known and the linear models apply only within the operating regions. 
Ibanez [5] takes a slightly different approach by assuming the system response to be 
periodic at the forcing frequency. An approximate transfer function is constructed. The 


276 


I 


tranfer function is dependent on the amplitude as well as on the exciting frequency and 
is valid only within the region of exciting frequencies. 

The direct equation approach is used by Yasuda, Kawamura and Watanabe [6], [7]. 
The input and output measurements of a dynamic process are expressed in a Fourier Series 
using, for example, an FFT algorithm. The system nonlinearity is represented as a sum of 
polynomials with unknown coefficients. Applying the principle of harmonic balance, die 
polynomial coefficients as well as the other system parameters are obtained accurately. 
Knowledge of the nonlinearity is needed to construct the polynomial. Truncation in the 
Fourier Series expansion of the input or output may lead to error. 

The regression estimation approach is used by Billings and Voon [8] and Greblick 
and Pawlak [9]. Billings and Voon [8] use the NARMAX model (Nonlinear Auto 
Regressive Moving Average model with exogenous inputs) to represent the nonlinear 
system. A stepwise regression method determines the significant terms in the NARMAX 
model. Then a prediction-error algorithm provides optimal estimates of the final model 
parameters. Greblick and Pawlak [9] represent the linear dynamic submodel by an ARMA 
model and the nonlineanties by a Borel function. A non-parametric kernel regression 
estimation is employed to obtain the final analytical model. 

Kortman and Unbehauen [10] and Distefano and Rath [11] use the minimization 
of an error cost function as a means of obtaining the coefficients of the functions used 
to represent the nonlinearities. The method presented by Kortman and Unbehauen [10] 
uses only system input and output information to estimate the polynomial representing 
the nonlinearities and the parameters of the linear components. It is robust in the 
presence of noise, although iteration is necessary. Distefano and Rath [11] present two 
techniques, a non-iterative direct identification and an iterative direct identification. In 
the first technique, measurement of all variables is required and the model parameters 
are obtained through the minimization of an error function. In the second technique, 
iteration is used to minimize a cost function yielding the system parameters in addition 
to the state trajectories. In Distefano and Rath [11], the nonlinear model form is also 
taken to be known. 


In other techniques, as in statistical linearization, a nonlinear relation is replaced by a 
linear equivalent gain. Broersen [12] extends the technique of statistical linearization by 
representing the nonlinearity as a linear combination of a number of arbitrary functions. 
Correlation techniques are then used to determine the coefficients of these functions. The 
number and type of functions selected depends on the desired accuracy as well as some 
knowledge of the system nonlinearity. Reasonable accuracy is obtained in the presence 
of noise and no iterations are necessary. Although some of the basic properties of the true 
nonlinear output are preserved, it is limited to only random excitation, and knowledge 
of all states and forcing terms is required. 

In the method of multiple scales (Hanagud, Mayyappa and Craig [13]), a perturbation 
solution to the nonlinear equation of motion is obtained. An objective function is built 
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employing an integral least squares approach. The minimization of the functional yields 
the unknown parameters. Data on only one field variable is necessary, and the method 
is effective in the presence of high noise. The method of multiple scales, however, is 
restricted to systems with small damping and slight nonlinearities and, as in most other 
methods, the form of the nonlinearity is assumed a priori. The method typically requires 
some algebraic manipulations which may be quite involved, and these manipulations are 
only valid for a particular assumed nonlinear form. If the assumed nonlinear form is 
changed, the algebra must be repeated. 

Another popular approach is to describe the nonlinear system using the Volterra or 
Wiener kernels. The Volterra series consists of the summation of impulse responses 
of increasing dimensionality. The Wiener series is also a set of orthogonal functions 
in which the input is white gaussian noise. Marmarelis and Udwadia [14], for example, 
estimate the first and higher order kernels appearing in the Volterra series using correlation 
techniques. Chen, Ishii and Suzumura [15] use cross-correlation functions in addition to 
the Volterra and Wiener series to describe nonlinear models and to show the relation 
between the system inner structure and the series. Although weakly nonlinear systems 
can be described by the first few kernels, for strongly nonlinear systems these series give 
accurate numerical results only at the expense of an excessive number of coefficients. 
This renders the analytical model impractical for control applications. 

Other popular series are orthogonal polynomials such as Legendre (Wang and Chan 
[16]), Chebyshev, and Jacobi (Horn and Chou [17]). Horn and Chou [17] expand the 
variables of the system into a shifted Jacobi series, reducing the nonlinear state equation 
into a linear algebraic matrix equation. The unknown parameters of the nonlinear system 
are then estimated using least squares. Even though the algorithm works well in the 
presence of noise, the nonlinear form must be known a priori. 

Another technique used for the identification of nonlinear systems is the extended 
Kalman filter. The extended Kalman filter is the linear Kalman filter applied to nonlinear 
systems by linearizing the nonlinear model into a Taylor series expansion about the 
estimated state vector. Yun and Shinozuka [18] apply the extended Kalman filter for the 
parameter estimation of a quadratic term. The state vector is augmented by including 
the unknown parameters in addition to the state variables. Through a series of iterations, 
the response, as well as the unknown parameters, are estimated by the Kalman filter. 
Among its disavantages are high sensitivity to initial conditions, in particular if the initial 
conditions are barely known. 

Hammond, Lo and Seager-Smith [19] use an optimal control technique based on 
optimal control methods employed for linear system deconvolution. The form of the 
linear model is assumed to be known as well as the input and the output. A cost functional 
consisting of the weighted sum of the square of the error (between the actual and estimated 
output) yields an optimal estimated input. The estimated input and the actual input are 
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used to obtain the nonlinearity as a function of the state variables. Although no previous 
assumption is made of the nonlinearities, there is no provision to deal with noise. 

In previous papers, the Minimum Model Error algorithm (MME) was explained 
in detail (Mook and Junkins [20]), modified for nonlinear identification (Mook [21]), 
and shown to accurately identify exotic nonlinearities in higher order systems (Stiy and 
Mook [22]). In this paper, it is shown how the MME algorithm successfully identifies 
nonlinearities using experimental data. An analytical model representing a harmonic 
oscillator with quadratic position feedback is studied. First, output data is obtained from 
a digital computer simulation of the nonlinear system and the quadratic term is identified 
to illustrate the accuracy of the technique on a known system. Second, an attempt is made 
to duplicate the nonlinear model using an analog computer. It is shown that despite the 
inability of the analog computer to produce a true quadratic term, the Minimum Model 
Error algorithm is capable of identifying a nonlinear model which accurately reproduces 
the analog output. The Minimum Model Error method produces a numerically stable 
identification regardless of the analog data initial conditions or record length. 


MME Algorithm 

In this section, we review the MME algorithm and how it is used to identify nonlinear 
dynamic systems. A more detailed explanation may be found in Mook and Junkins [20], 
Mook [21], and Stry and Mook [22]. 

The MME may be summarized as follows. Suppose there is a nonlinear system 
whose exact analytical representation is unknown, but for which output measurements 
are available. Using normal means (analysis, finite elements, etc.), a system model is 
constructed. As shown in [21]-[22], the MME will work well even if this system model 
is poor. The MME combines the assumed model with the measurements to determine the 
correct form of the nonlinear system. A correction term which represents the error in the 
model is added to the assumed model and a cost functional is formed. Minimization of 
the cost functional yields the model error. Subsequently, a least squares fit is performed 
on the error term to determine the correct form of the nonlinear system. 

Consider a forced nonlinear dynamic system which may be modeled in state-space 
form by the equation 

= m*) + m + /(*(<u(0) (i) 

where x(t) is the n x 1 state vector consisting of the system states, .4 is the n x n state 
matrix, F(t) is an n x 1 vector of known external excitation, and f(x(t),x(t)) is an 
n x 1 vector which includes all of the system nonlinearities. State -observable discrete 
time domain measurements are available for this system in the form 

y(ik) =£*(*(<*),<*) + £*, to<tk<tf (2) 
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where y(i k ) is an m x 1 measurement vector at time t k , g k is the accurate model of 
the measurement process, and v* represents measurement noise. v k is assumed to be 
a zero-mean, gaussian distributed process of known covariance R k . The measurement 
vector y(t k ) may contain one or more of the system states. To implement MME, assume 
that a model, which is generally not the true system model because of the difficulties 
inherent in obtaining the true system model, is constructed in state-vector form as 

x{t) = Ax(t) + F(t) (3) 

Here, we show a linear model because in practice, linearization is the most common 
approach to modeling nonlinear systems. MME uses the assumed linear model in (3) 
and the noisy measurements in (2) to find the model error. 

The model error, which includes the unknown nonlinear term(s) of the system, is 
represented by the addition of a correction term to the assumed linear model as 

x{t) = Ax(t) + F(t) + d(t) (4) 

where d(t) is the n x 1 correction term (dynamic model error) to be estimated later. 

A cost functional, J, that consists of the weighted integral square of the correction 
term plus the weighted sum square of the measurement-minus-estimated measurement 
residuals, is formed: 

J = - g k i^\ik)] T R^Wk) -£ fc (i(t*)><*)]} 

Jfc=l '■ 

+ T d{T) T Wd(r)dT ( 5 ) 

Jt o 

where M is the number of measurement times, x(t k ) is the estimated state vector and 
W is a weight matrix to be determined. 

j is minimized with respect to the correction term, d(t). The necessary conditions 
for the minimization lead to the following two point boundary value problem (TPBVP), 
(see Geering [23]), 


x(t) = Ax(t ) + jF(t) + d(t) 

(5a) 

A (t) = -A T \(t) 

(56) 

d(t) = 

(5c) 

\(t+) = \(t~) + 2iT Jb i£ 1 [y(4) -g k (m)^k)} 

(5i) 

Hk ^ it 

(5c) 

x(to) = Xo or X (to) = o 

x(tf) = xj or X{t f ) = 0 

(5/) 
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where A(f) is a vector of costates (Lagrange multipliers). Estimates of the states and 
of the dynamic model error are produced by the solution of this two-point boundary 
value problem. The estimates depend on the particular value of W. The solution is 
repeated until a value of W is obtained which produces state estimates which satisfy the 
“covariance constraint”, explained next. 

According to the covariance constraint, the measurement-minus-estimated measure- 
ment residual covariance matrix must match the measurement-minus-truth error covari- 
ance matrix. This may be written as 

(£(**) “&(£(**), <*)] T (y(<*) - 9 k (£(<*),<*)] ~ Rk (6) 

During the minimization, the weight W is varied until the state estimates satisfy the 
covariance constraint, i.e., the left hand side of Eq. (6) is approximately equal to 
the right hand side. The correction term or model error is, therefore, the minimum 
adjustment to the model required for the estimated states to predict the measurements 
with approximately the same covariance as the measurement error. 

The TPBVP represented by Eqs. (5a) to (5f) contains jumps in the costates and, 
consequently, in the correction term. As evident from Eq. (5d), the size of the jump is 
directly proportional to the measurement residual at each measurement time. The noisier 
the measurements, the larger the jump size. A multiple shooting algorithm, developed 
by Mook and Lew [24], converts this jump-discontinuous TPBVP into a set of linear 
algebraic equations which may be solved using any linear equation solver. Multiple 
shooting also facilitates the analysis of a large number of measurements, by processing 
the solution at the end of every set of jumps. 

After W has been determined such that the state estimates satisfy the covariance 
constraint, the final step in the identification procedure is to use a least squares algorithm 
to fit the model error d(t) to the unknown dynamic term(s). The error is expanded into 
some combination of linear and nonlinear terms, for example, 

d(t) = ax(t)+0x 2 (t) + yx 3 (t) + ... ( 7 ) 

where a, (3, 7, ... are unknown coefficients to be determined by least squares. Pre- 
sumably, the least squares fit of Eq. (7) will find zero coefficients for the terms in the 
expansion which are not part of the true model, and nonzero coefficients for the actual 
model terms. Eq. (7) may be sampled repeatedly to obtain 

d(ti) = ax(ti) + /?sr 2 (fi) + 7 * 3 (<i) + . . . 

<*(* 2 ) = ox(<2) + /3x 2 (t 2 ) + 7® 3 (<2) + . . . 

• • 

d(t\) = ax(ti) -f- px 2 (ti) + yx 3 (t t ) + ... 

or, in matrix form, 

£,xi = ixpILpxi (8) 
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where JP = [a fi 7 . . .] T is the vector of coefficients for the terms in d(t). Since 

estimates of d(t) are available continuously throughout the time domain, the parameter 
l may be chosen quite large to improve the least squares fit. Generally, because of the 
jump discontinuities in the model error estimates at the measurement times, it is desirable 
to pick the least squares sampling times in Eq. ( 8 ) at points other than the measurement 
times. The least squares estimate is found by minimizing the following cost functional 
with respect to P : 


$ = \D- MP) t [D - MP] (9) 

The solution is given by 

£ = (M t M)~ 1 M t D (10) 

The multiple shooting algorithm presented by Mook and Lew [24] was used to 
obtain the MME solutions used in the tests presented in this paper. It was assumed 
in the examples that MME obtained the dynamic error term without knowledge of the 
boundary conditions on £, so some distortion of the correction term at the initial and 
final times was expected due to the constraints of Eqs. (5e-5f), i.e., by assuming no state 
knowledge is available at to or if, we constrain A(<o) = 0 and A(ty) = 0 . Therefore, 
in all test cases, the initial and final ten percent of the correction term data was ignored 
in the least squares fit. 


Application Examples 

Two nonlinear equations of motion were studied, which represent the motion of 
an undamped harmonic oscillator with different amounts of quadratic position feedback 
(identical equations may arise in other physical systems as well). The equations in state 
space form are 


CM*, ;)0*u*.o 

(11) 


(12) 


where x is position, and the dot indicates differentiation with respect to time. No forcing 
was applied. 

In the following discussion, Eq. (11) is denoted Model A and Eq. (12) is denoted 
Model B. Different initial conditions were used for each system, for a total of five 
different tests. These are shown in Table 1. 
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Table 1. List of conditions used for each test 


Test# 

*(< o) 

*(<o) 

* 2 

Al 

0.00 

0.08 

-0.526 

A2 

0.00 

0.26 

-0.526 

A3 

0.00 

0.52 

-0.526 

B1 

0.00 

0.08 

-1.137 

B2 

0.00 

0.26 

-1.137 


To utilize MME, the linear part of Eqs. (11) and (12) was assumed to be known, 
rendering the model error equivalent to the nonlinear term, c * x 2 . Data for the 
MME nonlinear identification was generated from two different sources. First, noiseless 
position measurements were gathered from a digital computer simulation for all five tests. 
Application of MME to the measurements yielded an accurate estimate of the nonlinear 
term in each case. Then, an attempt was made to duplicate each system on an analog 
computer. Even though the analog computer did not reliably reproduce the quadratic 
term, the position measurements for all five tests were recorded and nonlinear models 
identified. MME proved capable of identifying accurate nonlinear models for the analog 
output. 

Digital computer simulation results 

One hundred noiseless position measurements were generated on a VAX 780 for the 
five test cases shown in Table 1.. A sampling rate of 4 Hertz was used. Three terms 
were employed in the least squares fit: x, x and x 2 . The resulting numerical values 
are shown in Table 2. 

Table 2. Least Squares estimates of the dynamic model error employing analyt- 
ically generated measurements in the MME algorithm. 


Test 

MME 

MME 

MME 

Analytic 

Analytic 

Analytic 

# 

X 

X 

x 2 

X 

X 

x 2 

Al 


4.90 xlO -4 

-0.526 

0.0 

0.0 

-0.526 

A2 


-3.00xl0 -4 

-0.527 

0.0 

0.0 

-0.526 

A3 



-0.528 

0.0 

0.0 

-0.526 

B1 



-1.138 

0.0 

0.0 

-1.137 

B2 

eeeeh 

-6.80xl0 -3 

-1.138 

0.0 

0.0 

-1.137 


MME identifies the quadratic term with great accuracy in all five tests. A plot of 
the estimated, analytical and measured position is shown in Figure (la) for test case Al. 
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ESTIMATED**). ANALYTICA L(-) ANO MEASURED POSITION 


Figure (lb) presents the predicted model error and the dynamic model error estimated 
by MME for test case Al. MME estimates the position and correction term with great 
accuracy. Similar results are obtained for tests A2, A3, Bl, and B2, but are not shown. 




(a) (b) 


Figure 1 (a) Analytical, measured (A), and MME estimated position 
for test Al using digital computer simulated measurements. The 
analytical and MME estimates are essentially identical (solid line), 
(b) MME estimate (+) and actual model error. 


Analog computer results 

Three hundred fifty position measurements were generated on a Comdyna GP-6 
analog computer for all five test cases. One hundred measurements with a sampling rate 
of 4 Hertz were used in the analysis. The identification procedure yielded the numerical 
values shown in Table 3. 
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Table 3. Least Square estimates of the nonlinear terms using measurements 
generated by the analog computer. 


Test 

# 

MME 

X 

MME 

X 

MME 

* 2 

A1 

-0.20 

1.04xl0 -2 

-6.17 

A2 

-5.06 xlO -2 

-1.89xl0 -3 

-1.322 

A3 

-6.44 xlO -3 

Q3BE9 

-0.689 

B1 

0.10 

EZHEH1 

-3.47 

B2 

2.55 xl(T 2 

7.42 xlO" 3 

-1.265 


The numerical results for the least squares lit of the error term did not match the 
analytically predicted coefficients. The reason for the numerical discrepancy was the 
analog s failure to produce a dependable quadratic term. Table 4 shows some position 
values squared by the analog. The analog consistently produced an error in the quadratic 
term. The recorded data, although containing errors due to quadratic term, is believed 
to be practically noiseless. 

Table 4. Quadratic term produced by the analog computer. 


X 

* 2 

Analog 

2.00 

4.00 

4.30 

2.50 

6.25 

7.00 

3.00 

9.00 

9.50 


Figures (2)- (6) show the analytical position, analog measurements and position 
predicted by the MME analysis for all analog tests. The MME identification produced 
good state estimates and a model which matched the measured data much better than the 
analytical models in Eqs. (11) and (12). 


Note, these results were obtained without knowledge of the initial or final state vector 
value. As shown in Eqs. (5e) and (5f), by setting the initial and final costate values to 
zero, no knowledge of the initial or end conditions are necessary. Also, the same results 
presented in Table 3 were obtained when using all three hundred and fifty measurements 
instead of one hundred measurements in the MME algorithm. 
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Figure 2 Analytical, measured (A), and MME estimated position 
for test A1 using analog computer measurements. The 
MME estimates are essentially identical to the measurements 



Figure 3 Analytical, measured (A), and MME estimated position 
for test A2 using analog computer measurements. The 
MME estimates are essentially identical to the measurements 




Figure 4 Analytical, measured (A), and MME estimated position 
for test A3 using analog computer measurements. The 
MME estimates are essentially identical to the measurements 



TlMt(S) 


Figure 5 Analytical, measured (A), and MME estimated position 
for test B1 using analog computer measurements. The 
MME estimates are essentially identical to the measurements 
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Figure 6 Analytical, measured (A), and MME estimated position 
for test B2 using analog computer measurements. The 
MME estimates are essentially identical to the measurements 

Conclusion 

In this paper, an MME — based algorithm was used to identify the quadratic term of a 
nonlinear harmonic oscillator. For demonstration purposes, data was obtained from two 
sources. Output data obtained from a digital computer simulation was used to verify the 
accuracy of the method. Then, data from an analog computer was used as a test of the 
method on “real” data. It was shown that despite the inability of the analog computer 
to reproduce a true quadratic term, the MME algorithm was capable of identifying a 
nonlinear model which accurately reproduced the analog output. This result indicates 
that the method is robust with respect to (lack of) a priori knowledge of the system 
dynamics. The identification was accurate regardless of initial conditions or data record 
length, indicating that the method is also robust with respect to those variables. 
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